/*******************************************************************************

This code file produces Table A12, "WTP and Cost by Quartile of MVPF Distribution."

programs to run:
-metafile.do
-run_program.ado
-est_life_impact_nyc.ado
-get_wtp_housing.ado

*******************************************************************************/

*** Manage settings

	run "$dir/code/modules/settings.do"
	
	do "$code/modules/stata-tex.do"
	cd "$tables/estimates"
	
********************************************************************************
* Tempfile to store results
********************************************************************************

	clear
	
	local NTAs = 179
	local B = 100
	local N = `B'*`NTAs'
	
	set obs `N'
	
	gen b = .
	gen nta = . 
	gen nta_code = ""
	gen NTAName = ""
	gen mvpf = .
	gen wtp = .
	gen mcost = .
	gen program_cost = .
	gen cost = .
	gen incidence = .
	
	save "$data/clean/mvpf_by_nta_boot.dta", replace
	
	forvalues b = 1/`B' {
	
	di `b'
	
	quietly {
	
********************************************************************************
* Obtain marginal cost per unit by NTA
********************************************************************************

	use "$data/clean/cleaned_data.dta", clear
	
	bsample, cluster(nta)
	tempfile boot
	save `boot', replace

	logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
	mat b = e(b)
	local coef_dtaxrate = b[1,1]
	
	predict xb, xb
	gen pr = exp(xb)/(1+exp(xb))
	gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
	
	gen mvalue_bldg = (assesstot - assessland)/(frac_assess*underassess)
		
	gen assessprop = pr * (assesstot - assessland) * dtaxrate_onsite/sh_taxable
	gen assessprop_plus = pr_plus * (assesstot - assessland) * (dtaxrate_onsite/sh_taxable + 0.001)
	
	local sh_inclusionary = 0.20
	gen inclusionaryunits = `sh_inclusionary' * pr * unitsres 
	gen inclusionaryunits_plus = `sh_inclusionary' * pr_plus * unitsres
	
	collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus (mean) xcoord_latlong ycoord_latlong, by(nta_code nta)
		
	gen mcostperunit = (1-$beta) * (assessprop_plus-assessprop) / (inclusionaryunits_plus-inclusionaryunits)
	
	tempfile mcost
	save `mcost', replace
	
	use "$data/GIS_boundaries/neighborhood_boundaries/df.dta", clear
	merge 1:1 nta_code using `mcost'
	
	keep nta NTAName nta_code mcost
	
	save `mcost', replace
	
********************************************************************************
* Obtain 421-a incidence by NTA
********************************************************************************
	
	use `boot', clear
	
	logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
		
	mat b = e(b)
	local coef_dtaxrate = b[1,1]
	predict xb, xb
	
	gen dlprofit = (1 / `coef_dtaxrate') * log(1 + exp(xb))
	gen mvalue = (assesstot-assessland)/(frac_assess*underassess)

	collapse (mean) dlprofit dtaxrate_onsite (rawsum) mvalue [aw=mval], by(borough nta nta_code)
	
	gen incidence = dlp/dtax
	
	replace incidence = 1 if incidence > 1 & !missing(incidence)
	
	bys borough: gegen inci_mean = mean(incidence) [aw=mvalue]
	replace incidence = inci_mean if missing(incidence)
	drop inci_mean borough
	
	keep nta_code incidence
	drop if missing(nta_code)
	
	tempfile incidence
	save `incidence', replace
	
********************************************************************************
* Merge together and run MVPF program
********************************************************************************
	
	preserve
	use "$data/raw/oi_pred_effect_by_nta_boot.dta", clear
	
	keep if b == `b'
	
	tempfile oi
	save `oi', replace
	restore
	
	merge 1:m nta_code using `oi', nogen
	merge m:1 nta_code using `mcost', nogen
	merge m:1 nta_code using `incidence', nogen
	
	keep if !missing(mcostperunit) & !missing(incidence) & !missing(pct_diff)
	
	local n = _N
	
	forvalues i = 1/`n' {
						
		global year_cost = mcost[`i']
		global pct_diff = pct_diff[`i']
		global incidence = incidence[`i']
		
		local ntaname = NTAName[`i']
		local nta = nta[`i']
		local nta_code = nta_code[`i']
		
		local mcost = mcost[`i']
		local incidence = incidence[`i']
	
		quietly run_program nyc421a
		local mvpf = $MVPF_nyc421a
		local wtp = $WTP_nyc421a
		local program_cost = $program_cost_nyc421a
		local cost = $cost_nyc421a
	
		preserve
		use "$data/clean/mvpf_by_nta_boot.dta", clear
	
		quietly replace mvpf = `mvpf' if _n == `i' + (`b'-1)*`NTAs'
		quietly replace wtp = `wtp' if _n == `i' + (`b'-1)*`NTAs'
		quietly replace program_cost = `program_cost' if _n == `i'+ (`b'-1)*`NTAs'
		quietly replace cost = `cost' if _n == `i' + (`b'-1)*`NTAs'
		
		quietly replace mcost = `mcost' if _n == `i' + (`b'-1)*`NTAs'
		quietly replace incidence = `incidence' if _n == `i' + (`b'-1)*`NTAs'
		
		quietly replace b = `b' if _n == `i' + (`b'-1)*`NTAs'
		quietly replace NTAName = "`ntaname'" if _n == `i' + (`b'-1)*`NTAs'
		quietly replace nta = `nta' if _n == `i' + (`b'-1)*`NTAs'
		quietly replace nta_code = "`nta_code'" if _n == `i'+ (`b'-1)*`NTAs'
		
		quietly save "$data/clean/mvpf_by_nta_boot.dta", replace
		
		restore
	
	}
	}
	}
	
	use "$data/clean/mvpf_by_nta_boot.dta", clear
	
	drop if missing(nta_code)
			
	replace mvpf = 0 if wtp < 0
	replace mvpf = 9999 if mvpf < 0 & wtp > 0 & !missing(wtp)
	
	format mvpf %6.2f
	
	save "$data/clean/mvpf_by_nta_boot.dta", replace
	
********************************************************************************
* Run MVPF program for quartile cuts
********************************************************************************
	
	use "$data/clean/mvpf_by_nta.dta", clear
	drop if missing(nta_code)
	
	gegen rank_mvpf = rank(mvpf)
	
	summ rank_mvpf
	local max = r(max)
	replace rank_mvpf = rank_mvpf / `max'
	
	keep nta_code rank_mvpf
	egen quartiles = cut(rank_mvpf), group(4)
	
	merge 1:m nta_code using "$data/clean/mvpf_by_nta_boot.dta", nogen
	
	collapse (mean) wtp_mean=wtp cost_mean=cost program_cost_mean=program_cost ///
			(semean) wtp_se=wtp cost_se=cost program_cost_se=program_cost, ///
			by(quartiles)
		
	forvalues q = 1/4 {
	
		local wtp_q`q'_est = wtp_mean[`q']	
		local wtp_q`q'_se = wtp_se[`q']
		
		local cost_q`q'_est = cost_mean[`q']
		local cost_q`q'_se = cost_se[`q']

		local gcost_q`q'_est = program_cost_mean[`q']
		local gcost_q`q'_se = program_cost_se[`q']
	
		insert_into_file using analysis_t21.csv, key(wtp_q`q'_est) value(`wtp_q`q'_est') format(%12.0fc)
		insert_into_file using analysis_t21.csv, key(wtp_q`q'_se) value(`wtp_q`q'_se') format(%12.0fc)
		
		insert_into_file using analysis_t21.csv, key(cost_q`q'_est) value(`cost_q`q'_est') format(%12.0fc)
		insert_into_file using analysis_t21.csv, key(cost_q`q'_se) value(`cost_q`q'_se') format(%12.0fc)
		
		insert_into_file using analysis_t21.csv, key(gcost_q`q'_est) value(`gcost_q`q'_est') format(%12.0fc)
		insert_into_file using analysis_t21.csv, key(gcost_q`q'_se) value(`gcost_q`q'_se') format(%12.0fc)
		
	}
	
********************************************************************************
* Create table
********************************************************************************
	
	cat analysis_t21.csv

	cap erase "$tables/output/analysis_t21.tex"

	cd "$code/modules"
	capture table_from_tpl, t("$tables/templates/analysis_t21.tex") ///
							r("$tables/estimates/analysis_t21.csv") ///
							o("$tables/output/analysis_t21.tex") 
			
	capture table_from_tpl, t("$tables/templates/analysis_t21.tex") ///
							r("$tables/estimates/analysis_t21.csv") ///
							o("$tables_overleaf/analysis_t21.tex") 
							
	exit
